6  Visualização de Dados Geográficos

Conteúdos

Visualização de dados resultantes de modelos de machine learning e de dados geográficos.

6.1 Dados Geográficos

Sistema de Informação Geográfica (SIG) é um sistema que permite manipular, armazenar, e fazer a análise espacial de dados geográficos.

6.1.1 Tipo de dados

pontos - Point(2,10) linha - LineString([(1,2),(1,5),…]) polígonos - Polygon([(13,1),(14,4),…])

Preparação do ambiente python:

através do prompt ANACONDA, instalar os packages: GEOPANDAS FOLIUM NBB

### Coordenadas Geográficas

Latitude e Longitude

Conhecer a Projeção Geográfica é importante para combinar dados de diferentes fontes de informação.

Coordinate Reference System (CRS): define como as coordenadas são representadas no plano.

EPSG Geodetic Parameter Dataset: registo global dos diferentes CRS.

EPSG: European Petroleum Survey Group para identificar os diferentes sistemas de coordenadas.

WGS84 - World Geodetic System 1984 associado ao sistema de posicionamento global (GPS).

• EPSG 4326
• coordenadaas em graus decimais

Pseudo-Mercator/WGS84

• EPSG 3857
• coordenadas em metros

6.1.1.1 Exercicio

Mostrar Influencia do tipo de Projeção Geográfica

# Este Script Permite Visualizar diferentes Projeções
import geopandas as gpd
import matplotlib.pyplot as plt

# Definir Figura do MATPLOTLIB
# Conseguem mostrar mais SUBPLOTS alterando argumeto 2 (1,2)
fig, ax = plt.subplots(1, 2, figsize=(15, 5))

# Carregar a geogrfia do mundo a partir d eum shapefile dentro um ficheiro ZIP
# Atenção Mudar o caminho!
world = gpd.read_file(r"data\ne_110m_admin_0_countries.zip")

# # Alternativa - Ler os Dados de um Link GEOJSON:
# url = "http://d2ad6b4ur7yvpq.cloudfront.net/naturalearth-3.3.0/ne_110m_land.geojson"
# # Ler Ficheiro GEOJSON do URL
# world = gpd.read_file(url)

# Sistema de coordenadas definidos na listagem em baixo:
# 
# WGS 84 (Lat/Lon): ['EPSG:4326','WGS 84 (Lat/Lon)']
# Web Mercator Metros: ['EPSG:3857','Web Mercator']
# Vander Grinten: ['ESRI:54029','World Robinson']
# Robinson: ['ESRI:54030','World Robinson - National Geographic']
# Gall-Peters - ['ESRI:53016','Gall Stereographic']
# Peirce Quincuncial: ['ESRI:54091','Peirce quincuncial North Pole in a square']
# winkel-Tripel: ['ESRI:53018','Sphere Winkel I']
# Goode Homolosine (Molweide) - ['EPSG:7619','Interrupted Goode Homolonsine']
# Albers Equal Area: ['ESRI:102008','North America Albers Equal Area Conic']
# Eckert: ['ESRI:53010','Sphere Eckert VI']
# Portugal Continental: ['EPSG:3763','Portugal Continental: PT-TM06/ETRS89']
# Açores UTM Fuso 26: ['EPSG:5015','TRF93 / PTRA08 - UTM fuso 26 - Grupo Central e Oriental do Arquipélago dos Açores']

# sistema de Coordenadas para Visualizar - 12 exemplos:
sistcoord = [['EPSG:4326','WGS 84 (Lat/Lon)'], ['EPSG:3857','Web Mercator'],
             ['ESRI:54029','World Robinson'], ['ESRI:54030','World Robinson - National Geographic'], 
             ['ESRI:53016','Gall Stereographic'], ['ESRI:54091','Peirce quincuncial North Pole in a square'],
             ['ESRI:53018','Sphere Winkel I'],['EPSG:7619','Interrupted Goode Homolonsine'],
             ['ESRI:102008','North America Albers Equal Area Conic'], ['ESRI:53010','Sphere Eckert VI'],
             ['EPSG:3763','Portugal Continental: PT-TM06/ETRS89'], ['EPSG:5015','TRF93 / PTRA08 - UTM fuso 26 - Grupo Central e Oriental do Arquipélago dos Açores']]


# Para mostrar
mapa1 = 0
mapa2 = 11

# ------------------------------------------
# Mudar Projeçoes na indicação da Listagem
#world.to_crs(epsg=2264).plot(ax=ax[0])

world.to_crs(sistcoord[mapa1][0]).plot(ax=ax[0])
ax[0].set_title(sistcoord[mapa1][1])

#world.to_crs(epsg=4087).plot(ax=ax[1])
world.to_crs(sistcoord[mapa2][0]).plot(ax=ax[1])
ax[1].set_title(sistcoord[mapa2][1])


plt.tight_layout()
plt.show()

6.1.2 GeoPandas e MatPlotLib

As mesmas funcionalidades do package Pandas podem ser usadas com o GeoPandas, usa também outros packages como Shapely e Fiona

6.1.3 Criar dados GeoPandas

  # Importar GeoPandas
import geopandas as gpd

# Carregar Dados com read_file (Shapefile preferivel, GeoJSON -  Mais Lento)
# Mudar Caminho para onde estão os dados - atenção de ter os ficheiros .shp\.shx\.dbf\.prj
file_path = r"data\NUTS3_2015_PT.shp"

# Definir o encoding para evitar problemas de desenho dos nomes
encoding = 'utf-8'  
# Ler Shapefile:
gdf_nuts3 = gpd.read_file(file_path, encoding=encoding)

print(gdf_nuts3.info())
print(gdf_nuts3.head())
gdf_nuts3.loc[1,'geometry']

#gdf_nuts3.head()

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 25 entries, 0 to 24
Data columns (total 8 columns):
 #   Column      Non-Null Count  Dtype   
---  ------      --------------  -----   
 0   NUTS1       25 non-null     object  
 1   NUTS2       25 non-null     object  
 2   NUTS3       25 non-null     object  
 3   NUTS3_DSG   25 non-null     object  
 4   SUM_AREA_K  25 non-null     float64 
 5   Shape_Leng  25 non-null     float64 
 6   Shape_Area  25 non-null     float64 
 7   geometry    25 non-null     geometry
dtypes: float64(3), geometry(1), object(4)
memory usage: 1.7+ KB
None
  NUTS1 NUTS2 NUTS3                    NUTS3_DSG   SUM_AREA_K     Shape_Leng  \
0     1    11   111                   Alto Minho  4005.236217  348704.286332   
1     1    11   112                       Cávado  2230.779298  308604.284583   
2     1    11   119                          Ave  2589.276604  418041.964330   
3     1    11   11A  Área Metropolitana do Porto  3596.793275  523879.094897   
4     1    11   11B                  Alto Tâmega  5240.318314  473818.956252   

     Shape_Area                                           geometry  
0  4.006183e+09  MULTIPOLYGON (((-911484.585 5182506.803, -9118...  
1  2.231519e+09  POLYGON ((-895962.526 5133087.503, -895581.858...  
2  2.589188e+09  POLYGON ((-893252.407 5115827.796, -892566.222...  
3  3.597098e+09  POLYGON ((-974333.659 5081342.266, -973571.650...  
4  5.239871e+09  POLYGON ((-880130.085 5150064.989, -878083.752...  

Mostrar Objecto Geometry

# Informação de um Polygon
# print(gdf_nuts3.loc[1,'geometry'])

Objecto de Geometria

# Tipo Objecto GeoPandas
print (type(gdf_nuts3))
# Tipo de Dados da coluna de Geometria
print (type(gdf_nuts3.geometry))

<class 'geopandas.geodataframe.GeoDataFrame'>
<class 'geopandas.geoseries.GeoSeries'>

Visualizar GDF

# Mostrar Geografia
#import geopandas as gpd
import matplotlib.pyplot as plt

gdf_nuts3.plot(column = 'NUTS3_DSG',
              legend = False)
plt.show()

Selecionar registos dum GDF

# Selecção utilizando os operadores comuns
# Selecção pode ser efetuada com ou sem a função LOC
gdf_nuts3_continente = gdf_nuts3[gdf_nuts3['NUTS3'] < '2']
gdf_nuts3_sel = gdf_nuts3.loc[gdf_nuts3.NUTS2 == '11']

# Mostrar numero registos selecionados
print("Nº registos Continente", len(gdf_nuts3_continente))
print("Nº registos NUTS2 11", len(gdf_nuts3_sel))

# Selecionar utilizando comparação de um String
selected_rows_like = gdf_nuts3[gdf_nuts3['NUTS3'].str.startswith('11')]
print('Nº registos que começam com 11',len(selected_rows_like))
                        
# Selecionar Linhas quando existem numa Listagem
lista_selecao = ['11A', '11E','16E','170']
selected_rows = gdf_nuts3[gdf_nuts3.NUTS3.isin(lista_selecao)] 
print('Nº registos que estão numa lista - vs 1',len(selected_rows))
selected_rows = gdf_nuts3.loc[gdf_nuts3.NUTS3.isin(lista_selecao)] 
print('Nº registos que estão numa lista - vs 2',len(selected_rows))

# Selecao dos Dados a Mostrar
gdf_nuts3_sel = gdf_nuts3.loc[gdf_nuts3.NUTS2 == '11']
gdf_nuts3_sel.plot(column = 'NUTS3_DSG',
              legend = True,
                  edgecolor = 'black')

Nº registos Continente 23
Nº registos NUTS2 11 8
Nº registos que começam com 11 8
Nº registos que estão numa lista - vs 1 4
Nº registos que estão numa lista - vs 2 4

<Axes: >

Visualização Interactiva com Explore

#import geopandas as gpd

# Selecao dos Dados a Mostrar
gdf_nuts3_sel = gdf_nuts3.loc[gdf_nuts3.NUTS2 == '11']
gdf_nuts3_sel.explore(column = 'NUTS3_DSG',
              legend = True,
                  edgecolor = 'black')

Make this Notebook Trusted to load map: File -> Trust Notebook

Mostrar geografia

# Mostrar Geografia

# Mostrar Áreas NUTS3 Região Centro
gdf_nuts3_sel = gdf_nuts3.loc[gdf_nuts3.NUTS2 == '16']

# Definir Legenda 
# Definir Localização (loc), Localização em relação ao eixos (bbox_to_anchor), Nº colunas (ncol)
lgnd_kwds = {'title': 'Nuts3 (Centro)',
               'loc': 'upper left', 
             'bbox_to_anchor': (1, 1.03), 
             'ncol': 1}
# Argumetnos: Atributo (column), Colormap (cmap), colocar Legenda (legend), 
#  Legenda definida (legend_kwds), 
gdf_nuts3_sel.plot(column = 'NUTS3_DSG',
              cmap = 'tab20',
              legend = True,
              legend_kwds  = lgnd_kwds,
              edgecolor = 'dimgray')

# Rotulos  Eixos
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Áreas NUTS3 de Portugal Continental')

plt.show()

6.1.3.0.1 Exercício

# packages necessários
import geopandas as gpd

# Carregar Dados com read_file (Shapefile preferivel, GeoJSON -  Mais Lento)
# Mudar Caminho para onde estão os dados - atenção de ter os ficheiros .shp\.shx\.dbf\.prj
file_path = r"data\CAOP20200_MN_PT.shp"

# Definir o encoding para evitar problemas de desenho dos nomes
encoding = 'utf-8'  
# Ler Shapefile:
gdf_caop = gpd.read_file(file_path, encoding=encoding)

print(gdf_caop.info())
print(gdf_caop.head())
print(gdf_caop.DISTRITO.unique()) 

gdf_caop.loc[1,'geometry']

print()
print(gdf_caop.NUTS3_02.unique()) # lista com a opções de NUTS3_02

# cria nova coluna para legenda
gdf_caop['DTMN2'] = gdf_caop['DTMNDSG'] + '(' + gdf_caop['DTMN'].astype(str) + ')'

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 308 entries, 0 to 307
Data columns (total 12 columns):
 #   Column      Non-Null Count  Dtype   
---  ------      --------------  -----   
 0   OBJECTID    308 non-null    int64   
 1   DTMN        308 non-null    object  
 2   DTMNDSG     308 non-null    object  
 3   ILHA        30 non-null     object  
 4   NUTS3_02    308 non-null    object  
 5   NUTS3_15    308 non-null    object  
 6   NUTS3_24    308 non-null    object  
 7   OBJECTID_2  308 non-null    int64   
 8   NUTS1_15    278 non-null    object  
 9   NUTS2_15    278 non-null    object  
 10  DISTRITO    278 non-null    object  
 11  geometry    308 non-null    geometry
dtypes: geometry(1), int64(2), object(9)
memory usage: 29.0+ KB
None
   OBJECTID  DTMN                DTMNDSG                ILHA NUTS3_02  \
0         1  4901                  Corvo       Ilha do Corvo      200   
1         2  4801       Lajes das Flores     Ilha das Flores      200   
2         3  4802  Santa Cruz das Flores     Ilha das Flores      200   
3         1  4201                  Lagoa  Ilha de São Miguel      200   
4         2  4202               Nordeste  Ilha de São Miguel      200   

  NUTS3_15 NUTS3_24  OBJECTID_2 NUTS1_15 NUTS2_15 DISTRITO  \
0    PT200      200           0     None     None     None   
1    PT200      200           0     None     None     None   
2    PT200      200           0     None     None     None   
3      200      200           0     None     None     None   
4      200      200           0     None     None     None   

                                            geometry  
0  POLYGON ((-3463610.566 4817991.232, -3463602.6...  
1  POLYGON ((-3479397.272 4791980.698, -3479270.4...  
2  POLYGON ((-3474231.695 4797050.671, -3474205.1...  
3  POLYGON ((-2845070.611 4547832.086, -2845103.7...  
4  POLYGON ((-2812174.924 4560067.448, -2812097.9...  
[None 'Évora' 'Lisboa' 'Portalegre' 'Porto' 'Santarém' 'Castelo Branco'
 'Faro' 'Viana do Castelo' 'Viseu' 'Coimbra' 'Setúbal' 'Vila Real'
 'Aveiro' 'Braga' 'Guarda' 'Bragança' 'Beja' 'Leiria']

['200' '183' '182' '16B' '115' '114' '185' '166' '150' '111' '165' '162'
 '16C' '181' '117' '118' '161' '116' '113' '164' '168' '167' '171' '172'
 '16A' '169' '184' '163' '112' '300']

vizualisar com explore

gdf_caop.explore()

Make this Notebook Trusted to load map: File -> Trust Notebook

contagens e filtros

# conta o número de registos por NUTS3_02
print(gdf_caop.NUTS3_24.value_counts())

# selecionar linhas a partir de uma lista
lista_selecao = ['11A','11E','16E','170']

select_rows = gdf_caop.loc[gdf_caop.NUTS3_24.isin(lista_selecao)]

print('Nº de registos que estão na lista ', lista_selecao, len(select_rows))

NUTS3_24
200    19
11D    19
192    19
11A    17
150    16
111    16
1C3    15
196    15
194    14
1C4    14
1C2    13
1D1    12
191    11
1D2    11
1D3    11
11C    11
300    11
193    10
11E     9
1A0     9
1B0     9
195     8
119     8
11B     6
1C1     5
Name: count, dtype: int64
Nº de registos que estão na lista  ['11A', '11E', '16E', '170'] 26

Fazer a Seleção de um DISTRITO e mostrar os munícipios

import matplotlib.pyplot as plt

# Selecao dos Dados a Mostrar
gdf_caop_sel = gdf_caop.loc[gdf_caop.DISTRITO == 'Guarda'] 

# Definir Legenda 
# Definir Localização (loc), Localização em relação ao eixos (bbox_to_anchor), Nº colunas (ncol)
lgnd_kwds = {'title': 'Municipios',
              'loc': 'upper left', 
             'bbox_to_anchor': (1, 1.03), 
             'ncol': 1}
             
gdf_caop_sel.plot(column = 'DTMN2',
              legend = True,
              legend_kwds  = lgnd_kwds,
              edgecolor = 'black')

# Rotulos  Eixos
plt.xlabel('X Coordinados')
plt.ylabel('Y Coordinados')
plt.title('Municipios do Distrito da Guarda')
# Rodar xticks; ha é o alinhamento
plt.xticks(rotation=45, ha='right')
              


# **Codigo Alternativo Utilizando ax objecto**
# fig, ax = plt.subplots()
# gdf_caop_sel.plot(ax=ax,
#               column = 'DTMN2',
#               cmap = 'tab20',
#               legend = True,
#               legend_kwds  = lgnd_kwds,
#               edgecolor = 'dimgray')
# 
# # Rotulos  Eixos - Diferentes Funções
# ax.set_xlabel('X Coordinados')
# ax.set_ylabel('Y Coordinados')
 
plt.show()

usando o explore

gdf_caop_sel.explore(column = 'DISTRITO',
              legend = True,
                  edgecolor = 'black')

Make this Notebook Trusted to load map: File -> Trust Notebook

6.1.4 Importar dados de tabelas

6.1.4.1 ACES

# Importar dados DGS Utentes em Centros de Saude
# Ver referencia:  https://dados.gov.pt/pt/datasets/evolucao-mensal-de-utentes-atendidos-nos-centros-de-saude-agregado-por-aces-no-ambito-da-saude-oral-nos-cuidados-de-saude-primarios-socsp/#resources
# Atributos: Período, ARS, ACES, Localização Geográfica, Sexo, Faixa Etária, Nº Utentes,ID 
import pandas as pd

# Link DGS (Utentes Centro Saude)  
ficheiro = r'http://dados.gov.pt/pt/datasets/r/dc54ea6f-31f3-483b-a719-718d0d7451f3'
# Ler ficheiro do computador:
#ficheiro = r'C:\TEMp\utentes-atendidos-nos-centros-de-saude-no-ambito-da-soep.csv'

# Importar em DataFrame
encoding = 'utf-8'
df_utentes = pd.read_csv(ficheiro, sep=';', encoding=encoding)


# Mostrar informação df
print(df_utentes.head(5))
print(df_utentes.info())
print(df_utentes.describe())

   Período    ARS                           ACES Localização Geográfica  \
0  2019-07  Norte                   ULS Nordeste  41.8069684,-6.7587977   
1  2019-07  Norte  Douro I - Marão e Douro Norte  41.2968711,-7.7483727   
2  2019-07  Norte  Douro I - Marão e Douro Norte  41.2968711,-7.7483727   
3  2019-07  Norte           Douro II - Douro Sul  41.0953745,-7.8123805   
4  2019-07  Norte           Douro II - Douro Sul  41.0953745,-7.8123805   

        Sexo Faixa Etária  Nº Utentes  \
0   Feminino        20-34          30   
1  Masculino        35-49          23   
2   Feminino        50-64          70   
3  Masculino        20-34          41   
4  Masculino        35-49          80   

                                                  ID  
0                 2019-7/20-34/Feminino/ULS Nordeste  
1  2019-7/35-49/Masculino/Douro I - Marão e Douro...  
2  2019-7/50-64/Feminino/Douro I - Marão e Douro ...  
3        2019-7/20-34/Masculino/Douro II - Douro Sul  
4        2019-7/35-49/Masculino/Douro II - Douro Sul  
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 30692 entries, 0 to 30691
Data columns (total 8 columns):
 #   Column                  Non-Null Count  Dtype 
---  ------                  --------------  ----- 
 0   Período                 30692 non-null  object
 1   ARS                     30692 non-null  object
 2   ACES                    30692 non-null  object
 3   Localização Geográfica  30096 non-null  object
 4   Sexo                    30692 non-null  object
 5   Faixa Etária            30692 non-null  object
 6   Nº Utentes              30692 non-null  int64 
 7   ID                      30692 non-null  object
dtypes: int64(1), object(7)
memory usage: 1.9+ MB
None
         Nº Utentes
count  30692.000000
mean      59.380457
std       71.069367
min        1.000000
25%       15.000000
50%       36.000000
75%       78.000000
max      851.000000

passar para o GDF

# Verificar os Dados
# Total numero de Pontos Unicos para Ficheiro:
print('Num localizações Geográficos:',len(df_utentes['Localização Geográfica'].unique()))

# Verificar valores únicos de outras variáveis
print('Num ACES:',len(df_utentes['ACES'].unique()))
print('Num ARS:',len(df_utentes['ARS'].unique()))
print(df_utentes['Faixa Etária'].unique())
print(df_utentes['Sexo'].unique())

Num localizações Geográficos: 51
Num ACES: 79
Num ARS: 35
['20-34' '35-49' '50-64' '<20' '65 e +']
['Feminino' 'Masculino']

# DataFrame com dados Utentes: df_utentes

# Criar Colunas lat e Long a partir da coluna 'Localização Geográfica'
# no ficheiro original a informação existe em apenas uma coluna
print ("Tipo de atributo", type(df_utentes['Localização Geográfica']))
# Fazer um Split dos Valores utilizando virgula como separador
df_utentes[['lat', 'long']] = df_utentes['Localização Geográfica'].str.split(',', expand=True)

# Mostrar informação df
print(df_utentes.info())

Tipo de atributo <class 'pandas.core.series.Series'>
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 30692 entries, 0 to 30691
Data columns (total 10 columns):
 #   Column                  Non-Null Count  Dtype 
---  ------                  --------------  ----- 
 0   Período                 30692 non-null  object
 1   ARS                     30692 non-null  object
 2   ACES                    30692 non-null  object
 3   Localização Geográfica  30096 non-null  object
 4   Sexo                    30692 non-null  object
 5   Faixa Etária            30692 non-null  object
 6   Nº Utentes              30692 non-null  int64 
 7   ID                      30692 non-null  object
 8   lat                     30096 non-null  object
 9   long                    30096 non-null  object
dtypes: int64(1), object(9)
memory usage: 2.3+ MB
None

# Criar DataFrame com Pontos Unicos e Numero de Utentos

# Total numero de Pontos Unicos para Ficheiro:
print('Num loc:',len(df_utentes['Localização Geográfica'].unique()))

# Converter variaveis numericos para numeric
df_utentes['Nº Utentes'] = pd.to_numeric(df_utentes['Nº Utentes'], errors='coerce')
df_utentes['lat'] = pd.to_numeric(df_utentes['lat'], errors='coerce')
df_utentes['long'] = pd.to_numeric(df_utentes['long'], errors='coerce')

# Criar Novo DataFrame com ACES e Soma Nº Utentes
# Group by 'ARS', 'ACES', 'Localização Geográfica'
# O nº de utentes é um exemplo de um atributo que podemos utilizar para visualização
df_aces = df_utentes.groupby(['ARS', 'ACES', 'Localização Geográfica','lat', 'long']).agg({
    'Nº Utentes': 'sum'
}).reset_index()

print(df_aces.head())

Num loc: 51

        ARS                     ACES          Localização Geográfica  \
0  Alentejo         Alentejo Central           38.8442031,-7.5826619   
1   Algarve      Algarve I - Central  37.0274264,-7.9395983999999995   
2   Algarve      Algarve I - Central           37.0274264,-7.9395984   
3   Algarve  Algarve II - Barlavento           37.1387554,-8.5445093   
4   Algarve   AlgarveIII - Sotavento            37.383008,-7.7293275   

         lat      long  Nº Utentes  
0  38.844203 -7.582662       45349  
1  37.027426 -7.939598        5051  
2  37.027426 -7.939598       37627  
3  37.138755 -8.544509       35960  
4  37.383008 -7.729328       23136  

# Mesmo com Pandas DataFrame já é possível mostrar a geografia
# Importar Bibliotecas
#import pandas as pd
import matplotlib.pyplot as plt

# Mostrar a localização
plt.scatter(df_aces.long, df_aces.lat)
 
# Show the plot
plt.show()

# Criar GeoDataFrame 
#import pandas as pd
import geopandas as gpd
from shapely.geometry import Point

# Criar a Coluna de geometry no DF a partir informação atributos long e lat 
# Utiliza-se o Point Object do modulo shapely
# lambda - função anônima no Python 
# df_aces.apply - applica a função a cada linha do DF
# Argumento apply: axis = 1. Função deve ocorrer para todos as linhas
df_aces['geometry'] = df_aces.apply(lambda x: Point(float(x.long), float(x.lat)), axis=1)
 
# Create a GeoDataFrame from art and verify the type
gdf_aces = gpd.GeoDataFrame(df_aces, crs = "EPSG:4326", geometry = df_aces.geometry)

print(type(gdf_aces))

<class 'geopandas.geodataframe.GeoDataFrame'>

importar dados poligonos

# Shapefile NUTS3: NUTS3_2015_PT.shp (dados estão em WebMercator)
# Mostrar Geografia
#import geopandas as gpd
#import matplotlib.pyplot as plt

# Carregar Dados com read_file (Shapefile preferivel, GeoJSON -  Mais Lento)
# Mudar Caminho para onde estão os dados - atenção de ter os ficheiros .shp\.shx\.dbf\.prj
file_path = r"data\NUTS3_2015_PT.shp"

# Definir o encoding para evitar problemas de desenho dos nomes
encoding = 'utf-8'  
# Ler Shapefile:
gdf_nuts3 = gpd.read_file(file_path, encoding=encoding)

print(gdf_nuts3.info())

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 25 entries, 0 to 24
Data columns (total 8 columns):
 #   Column      Non-Null Count  Dtype   
---  ------      --------------  -----   
 0   NUTS1       25 non-null     object  
 1   NUTS2       25 non-null     object  
 2   NUTS3       25 non-null     object  
 3   NUTS3_DSG   25 non-null     object  
 4   SUM_AREA_K  25 non-null     float64 
 5   Shape_Leng  25 non-null     float64 
 6   Shape_Area  25 non-null     float64 
 7   geometry    25 non-null     geometry
dtypes: float64(3), geometry(1), object(4)
memory usage: 1.7+ KB

None

mostrar pontos com Polygons

# Existem diferenças entre o CRS!
print("CRS of the GeoDataFrame:", gdf_aces.crs)
print("CRS of the GeoDataFrame:", gdf_nuts3.crs)

# Diferença nos CRS -> Converter para WebMercator
gdf_aces_m = gdf_aces.to_crs("epsg:3857")

# Obter apenas Municipios do Continente:
gdf_nuts3_sel = gdf_nuts3.loc[gdf_nuts3.NUTS1 == '1']

# Mostrar os 2 mapas:
# Argumento cmap tem a predefinição de cores
# Lista Referencias ColorMaps: https://matplotlib.org/stable/gallery/color/colormap_reference.html
ax = gdf_nuts3_sel.plot(column='NUTS3', 
                        cmap='Set2', 
                        legend=True, 
                        figsize=(10, 8))

# Imprimir os pontos no mesmo ax (subplot)
# Lista MatPlotLib cores: https://matplotlib.org/stable/gallery/color/named_colors.html
gdf_aces_m.plot(ax=ax, 
                color='red',
                markersize=30, 
                edgecolor='Black', 
                label='Cidades')

# Adicionar Legenda
ax.legend()
plt.title('Localização das ACES por NUTS3')

# Mostrar o plot
plt.show()

CRS of the GeoDataFrame: EPSG:4326
CRS of the GeoDataFrame: EPSG:3857

O package Contextily permite adicionar um basemap

# Import contextily
import contextily
    
# Existem diferenças entre o CRS!
print("CRS of the GeoDataFrame:", gdf_aces.crs)
print("CRS of the GeoDataFrame:", gdf_nuts3.crs)

# Diferença nos CRS -> Converter para WebMercator
gdf_aces_m = gdf_aces.to_crs("epsg:3857")

# Obter apenas Municipios do Continente:
gdf_nuts3_sel = gdf_nuts3.loc[gdf_nuts3.NUTS1 == '1']

# Mostrar os 2 mapas:
# Argumento cmap tem a predefinição de cores
# Lista Referencias ColorMaps: https://matplotlib.org/stable/gallery/color/colormap_reference.html
# Adicionar Transperência (alpha) ao layer
ax = gdf_nuts3_sel.plot(column='NUTS3', 
                        cmap='Set2', 
                        legend=True, 
                        figsize=(10, 8),
                       alpha = 0.5)

# Imprimir os pontos no mesmo ax (subplot)
# Lista MatPlotLib cores: https://matplotlib.org/stable/gallery/color/named_colors.html
gdf_aces_m.plot(ax=ax, 
                color='red',
                markersize=30, 
                edgecolor='Black', 
                label='ACES')

# Adicionar um BaseMap no Ax
contextily.add_basemap(ax)

# Adicionar Legenda
ax.legend()
plt.title('Localização das ACES por NUTS3')

# Show the plot
plt.show()

CRS of the GeoDataFrame: EPSG:4326
CRS of the GeoDataFrame: EPSG:3857

Obter Informações sobre o Objecto de Geometria

A partir do geometry conseguimos obter um conjunto de caracteristicas

# Tipo de Dados da coluna de Geometria
print (type(gdf_aces.geometry))
print (type(gdf_nuts3_sel.geometry))
print(gdf_nuts3_sel.info())

<class 'geopandas.geoseries.GeoSeries'>
<class 'geopandas.geoseries.GeoSeries'>
<class 'geopandas.geodataframe.GeoDataFrame'>
Index: 23 entries, 0 to 22
Data columns (total 8 columns):
 #   Column      Non-Null Count  Dtype   
---  ------      --------------  -----   
 0   NUTS1       23 non-null     object  
 1   NUTS2       23 non-null     object  
 2   NUTS3       23 non-null     object  
 3   NUTS3_DSG   23 non-null     object  
 4   SUM_AREA_K  23 non-null     float64 
 5   Shape_Leng  23 non-null     float64 
 6   Shape_Area  23 non-null     float64 
 7   geometry    23 non-null     geometry
dtypes: float64(3), geometry(1), object(4)
memory usage: 1.6+ KB
None

obter Área e Centroid

# Obter Informaçao Geografia - Utilizando gdf_caop
# Importar geometry
from shapely.geometry import Point, Polygon

# Adicionar Nova coluna com area em Km2
# Atenção Deveriamos adicionar novos atributos sempre a (Geo)DataFrame de Origem
# Não adicionar a seleção, neste caso gdf_nuts3_sel!
gdf_nuts3['area_km2'] = gdf_nuts3['geometry'].area / 10**6  

# Ordenar por Area
area_gdf = gdf_nuts3.sort_values(by='area_km2', ascending=False)

# Print the GeoDataFrame with the new column
print(area_gdf[['NUTS3','NUTS3_DSG', 'area_km2']])

   NUTS3                     NUTS3_DSG      area_km2
19   184                Baixo Alentejo  13733.680581
22   187              Alentejo Central  12124.085355
16   16J     Beiras e Serra da Estrela  10917.782234
21   186                 Alto Alentejo  10140.463336
7    11E      Terras de Trás-os-Montes   9909.588021
18   181              Alentejo Litoral   8565.542590
8    150                       Algarve   7899.809228
14   16H                   Beira Baixa   7847.568545
11   16E             Região de Coimbra   7444.322853
6    11D                         Douro   7118.775462
20   185               Lezíria do Tejo   7110.526353
15   16I                    Médio Tejo   5639.532885
13   16G              Viseu Dão Lafões   5638.896615
4    11B                   Alto Tâmega   5239.870518
17   170  Área Metropolitana de Lisboa   4963.947294
12   16F              Região de Leiria   4157.670418
0    111                    Alto Minho   4006.182814
23   200    Região Autónoma dos Açores   3789.073859
9    16B                         Oeste   3711.889332
3    11A   Área Metropolitana do Porto   3597.098437
5    11C                Tâmega e Sousa   3236.885188
10   16D              Região de Aveiro   2942.595609
2    119                           Ave   2589.188487
1    112                        Cávado   2231.518979
24   300    Região Autónoma da Madeira   1136.433784

adicionar nova coluna com o centroide

# Criar Nova coluna com Centroid
# Adicionar Nova coluna com centroid
gdf_nuts3['mn_center'] = gdf_nuts3['geometry'].centroid  

# Print the GeoDataFrame with the new column
print(gdf_nuts3[['NUTS3','NUTS3_DSG', 'mn_center']])

   NUTS3                     NUTS3_DSG                         mn_center
0    111                    Alto Minho   POINT (-946980.726 5142757.347)
1    112                        Cávado   POINT (-941320.248 5103934.014)
2    119                           Ave   POINT (-909301.972 5085897.504)
3    11A   Área Metropolitana do Porto   POINT (-943820.268 5025611.477)
4    11B                   Alto Tâmega   POINT (-848339.419 5111654.442)
5    11C                Tâmega e Sousa   POINT (-905197.685 5040987.781)
6    11D                         Douro   POINT (-825322.310 5034791.692)
7    11E      Terras de Trás-os-Montes   POINT (-759535.812 5097352.556)
8    150                       Algarve   POINT (-905203.263 4473163.600)
9    16B                         Oeste  POINT (-1015558.750 4762587.544)
10   16D              Região de Aveiro   POINT (-949322.164 4958172.134)
11   16E             Região de Coimbra   POINT (-927707.382 4898122.619)
12   16F              Região de Leiria   POINT (-962098.688 4841198.729)
13   16G              Viseu Dão Lafões   POINT (-881906.240 4968423.075)
14   16H                   Beira Baixa   POINT (-825652.276 4850783.449)
15   16I                    Médio Tejo   POINT (-920812.797 4807647.369)
16   16J     Beiras e Serra da Estrela   POINT (-810641.235 4941054.650)
17   170  Área Metropolitana de Lisboa  POINT (-1006290.656 4683604.593)
18   181              Alentejo Litoral   POINT (-952473.033 4579022.311)
19   184                Baixo Alentejo   POINT (-870119.976 4561859.007)
20   185               Lezíria do Tejo   POINT (-958500.668 4738718.638)
21   186                 Alto Alentejo   POINT (-848325.426 4749055.249)
22   187              Alentejo Central   POINT (-872968.675 4665202.808)
23   200    Região Autónoma dos Açores  POINT (-3038400.829 4629675.418)
24   300    Região Autónoma da Madeira  POINT (-1887053.937 3862343.766)

calcular a distancia entre dois pontos

#import geopandas as gpd
#from shapely.geometry import Point, Polygon

# Calcular diancia entre 2 ACES ('Amadora' and 'Lisboa Central')

# Como dados são em Graus Decimais será necessário de os converter para Metricos
print("CRS gdf_aces:", gdf_aces.crs)

origem = 'Amadora'
destino = 'Grande Porto V - Porto Ocidental'
# Converter para WebMercator
gdf_aces_m = gdf_aces.to_crs("epsg:3857")

# Filtrar as geometrias para obter a Origem e Destino definidos
origem_geometry = gdf_aces_m[gdf_aces_m['ACES'] == origem]['geometry'].iloc[0]
destino_geometry = gdf_aces_m[gdf_aces_m['ACES'] == destino]['geometry'].iloc[0]

# Calculate the distance between 'Amadora' and 'Lisboa Central'
distance_km = origem_geometry.distance(destino_geometry) / 1000  # Converter para quilómetro

print(f"A distância entre {origem} e {destino} é aproximadamente  {distance_km:.2f} quilómetros.")

CRS gdf_aces: EPSG:4326
A distância entre Amadora e Grande Porto V - Porto Ocidental é aproximadamente  356.68 quilómetros.

6.1.5 Juntar Dados a um GDF

6.1.5.1 Censos 2021

importa dados de Censos2021 por NUTS3 e Municipio

# Obter Password e Utilizador para Ligacao SQL
#from getpass import getpass # para ler a password sem a mostrar
my_user = '"BRUNO.LIMA"'  
my_password = 'Pa$$w0rd5' # getpass("Password: ")

# Ler Dados da BD
# criar conexão
import cx_Oracle 
import pandas as pd
host = 'c21oradev01.int.ine.pt'
port = '1521'
service = 'FORMACAO'
dsn_tns = cx_Oracle.makedsn(host, port, service_name=service) 

# Criar a conexão com todos os elementos,
# incluingo user e password
conn = cx_Oracle.connect(user=my_user, password=my_password, dsn=dsn_tns) 

# Cursor:
# Criar o cursor na conexão conn que criámos antes
c = conn.cursor()

# ---------------------------------
# Ler View com Dados por NUTS3:

# SQL String
my_sql = """
select *
from BDIFRM.V_BGRI2021_N3_PT 
"""

# Executar o cursor c com a string como parâmetro
c.execute(my_sql)
# Criar Nomes colunas
names = [ x[0] for x in c.description]
df_n3_c2021 = pd.DataFrame(c.fetchall(), columns = names)

# Fechar cursor
c.close()

# fechar conexão
conn.close()

fazer JOIN dos dados GDF com as NUTS3

# Juntar Dados NUTS3 do C2021 a GDF com as areas NUTS3

# gdf_nuts3 obtido de Shapefile NUTS3: NUTS3_2015_PT.shp (dados estão em WebMercator)
# Colunas Comuns: NUTS3 e NUTS3 - podemos ver os valores em ambos os DF
#print(sorted(df_n3_c2021.NUTS3.unique()))
#print(sorted(gdf_nuts3.NUTS3.unique()))

# Fazer o Join, especificar: DF
gdf_nuts3_2 = gdf_nuts3.merge(df_n3_c2021, left_on='NUTS3', right_on='NUTS3', how='left') 

print(gdf_nuts3_2.info())

6.1.6 Criar Mapas Temáticos

Visualizar a classificação da relação entre 2 variáveis por NUTS3

# Import packages
# import matplotlib.pyplot as plt
# import pandas as pd
# import geopandas as gpd
 
# Exemplo criar Plot dos Areas NUTS3 - Mostrando total de população 65+ no total de população
# Utilizar geodataframe com join "gdf_nuts3_2"

# Selecionar Dados Portugal Continental:
# Fazer Seleção da NUTS1_x , houve rename da coluna apos join
gdf_nuts3_sel = gdf_nuts3_2.loc[gdf_nuts3_2.NUTS1_x == '1']

# Definir Legenda 
lgnd_kwds = {'loc': 'upper left', 
             'bbox_to_anchor': (1, 1.03), 
             'ncol': 2}

# Generate the choropleth and store the axis
# natural_breaks
ax = gdf_nuts3_sel.plot(column=gdf_nuts3_sel.N_INDIVIDUOS_65_OU_MAIS/gdf_nuts3_sel.N_INDIVIDUOS, 
                      scheme='quantiles', # natural_breaks, quantiles, equal_interval 
                      k=7, 
                      cmap='YlGn', 
                      legend=True,
                      edgecolor = 'dimgray',
                      legend_kwds  = lgnd_kwds)
 
# Remover frames, ticks e tick labels do axis
ax.set_axis_off()

plt.title('Relacao população 65+ no total de população')
plt.show()

6.1.6.1 Exercício

Criar Mapa Temático das NUTS com apresentação de cidades Estatisticas

# Fazer o Join de NUTS e dados dos Censos , especificar: DF
gdf_nuts3_2 = gdf_nuts3.merge(df_n3_c2021, left_on='NUTS3', right_on='NUTS3', how='left')

print(gdf_nuts3_2.info())

# Importar os dados das cidades
ficheiro = r'data\CIDADES_PONTOS_CONT.csv'

# Importar em DataFrame
encoding = 'utf-8'
df_cidades = pd.read_csv(ficheiro, sep=';', encoding=encoding)

df_cidades.head()

#from shapely.geometry import Point, Polygon

# Criar coluna com geometry (mudar nome df e nomes atributos)
df_cidades['geometry'] = df_cidades.apply(lambda x: Point(float(x.LONGITUDE), float(x.LATITUDE)), axis=1)

# Criar uma gdf a partir da coluna geometry (mudar nome df)
gdf_cidades = gpd.GeoDataFrame(df_cidades, crs = "EPSG:4326", geometry = df_cidades.geometry)

# Converter Cidades para WebMercator
gdf_cidades_m = gdf_cidades.to_crs("epsg:3857")

# Utilizar: gdf_nuts3_2
# Input variáveis: gdf_aces e gdf_nuts3_2

# Selecionar Dados Portugal Continental:
# Fazer Seleção da NUTS1_x , houve rename da coluna apos join
gdf_nuts3_sel = gdf_nuts3_2.loc[gdf_nuts3_2.NUTS1_x == '1']

# Definir Legenda 
lgnd_kwds = {'loc': 'upper left', 
             'bbox_to_anchor': (1, 1.03), 
             'ncol': 2}

# Generate the choropleth and store the axis
# natural_breaks
ax = gdf_nuts3_sel.plot(column=(gdf_nuts3_sel.N_EDIFICIOS_CONSTR_ANTES_1945+gdf_nuts3_sel.N_EDIFICIOS_CONSTR_1946_1980)/gdf_nuts3_sel.N_EDIFICIOS_CLASSICOS, 
                      scheme='quantiles', # natural_breaks, quantiles, equal_interval 
                      k=7, 
                      cmap='YlGnBu', 
                      legend=True,
                      edgecolor = 'gray',
                      legend_kwds  = lgnd_kwds)

# Adicionar os Pontos no mesmo AX
gdf_cidades_m.plot(ax=ax, color='lightcoral', markersize=10, edgecolor='Black', label='Cidades')

# Adicionar Legenda
plt.title('Rácio Edificios construidos até 1980 (com representação das cidades)')

# Remover frames, ticks e tick labels do axis
ax.set_axis_off()

plt.title('Relacao edificios até 1980 no total de edificios')
plt.show()

gdf_cidades_m.explore()

6.1.7 Folium

import folium

# Mostrar Zona do INE 
# Lisboa: 38.738345820101536, -9.138601637922605
# Porto: 41.150887961411634, -8.629046747840249
# Longitude e Latitude
ine = folium.Map(location = [41.150887961411634, -8.629046747840249], 
                 zoom_start = 15)

ine

Importar Dataset com NUTS3

# Importar GeoPandas
#import geopandas as gpd

# Carregar Dados com read_file (Shapefile preferivel, GeoJSON -  Mais Lento)
# Mudar Caminho para onde estão os dados - atenção de ter os ficheiros .shp\.shx\.dbf\.prj
file_path = r"data\NUTS3_2015_PT.shp"

# Definir o encoding para evitar problemas de desenho dos nomes
encoding = 'utf-8'  
# Ler Shapefile:
gdf_nuts3 = gpd.read_file(file_path, encoding=encoding)

gdf_nuts3.loc[1,'geometry']

Mostrar Localização Central da GDF das NUTS3

# Obter a Localização Central
# Print the head of the urban_polygon
#import geopandas as gdp
import folium

# Obter o centroid de toda a geometria
# Converter para 4326 e a seguir obter o centroid de união de toda a geometria
gdf_nuts3_sel = gdf_nuts3.loc[gdf_nuts3.NUTS1 == '1'] # Selecionar apenas o continente
centroid = gdf_nuts3_sel.to_crs(epsg=4326).unary_union.centroid

# Criar Listagem com localização de latitude longitude
center_map = [centroid.y, centroid.x]
# Criar Mapa e mostrar
folium_map = folium.Map(location=center_map, zoom_start=8, tiles='OpenStreetMap')

folium_map

Make this Notebook Trusted to load map: File -> Trust Notebook

adicionar geografia (GDF) ao mapa

# Adicionar a Geografia a mapa
# Print the head of the urban_polygon
# import geopandas as gdp
# import folium
# from shapely.geometry import Point, Polygon

# Converter para 4326 e a seguir obter o centroid de união de toda a geometria
centroid = gdf_nuts3_sel.to_crs(epsg=4326).unary_union.centroid

# Criar Listagem com localização de latitude  longitude
center_map = [centroid.y, centroid.x]
# Criar Mapa e mostrar
folium_map = folium.Map(location=center_map, zoom_start=7, tiles='OpenStreetMap')

# Adicionar Geografia folium map
# folium.GeoJson constructor
folium.GeoJson(gdf_nuts3_sel).add_to(folium_map)

folium_map

Make this Notebook Trusted to load map: File -> Trust Notebook

adicionar pop e tooltip

gdf_nuts3_sel.info()

# Incluir Informação de POPUP
# Print the head of the urban_polygon
#import geopandas as gdp
#import folium
import folium.plugins 

# Converter para 4326 e a seguir obter o centroid de união de toda a geometria
centroid = gdf_nuts3_sel.to_crs(epsg=4326).unary_union.centroid

# Criar Listagem com localização de latitude  longitude
center_map = [centroid.y, centroid.x]

# Criar Mapa e mostrar
folium_map = folium.Map(location=center_map, zoom_start=7, tiles='OpenStreetMap')

# Codigo Especifico para tooltip
tooltip = folium.GeoJsonTooltip(
    fields=["NUTS3"],
    aliases=["NUTS3"],
    localize=True,
    sticky=False,
    labels=True,
    style="""
        background-color: #F0EFEF;
        border: 2px solid black;
        border-radius: 3px;
        box-shadow: 3px;
    """,
    max_width=800,
)

# Definir popup e nome
folium.GeoJson(gdf_nuts3_sel,
        zoom_on_click = True,
        popup = folium.GeoJsonPopup(
            fields=['NUTS3','NUTS3_DSG'],
            aliases=['NUTS3','NUTS3_DSG']
        ),
        tooltip = tooltip       
               
    ).add_to(folium_map) 


folium_map
# folium_map.save(r'data\xxmap.html')

<class 'geopandas.geodataframe.GeoDataFrame'>
Index: 23 entries, 0 to 22
Data columns (total 8 columns):
 #   Column      Non-Null Count  Dtype   
---  ------      --------------  -----   
 0   NUTS1       23 non-null     object  
 1   NUTS2       23 non-null     object  
 2   NUTS3       23 non-null     object  
 3   NUTS3_DSG   23 non-null     object  
 4   SUM_AREA_K  23 non-null     float64 
 5   Shape_Leng  23 non-null     float64 
 6   Shape_Area  23 non-null     float64 
 7   geometry    23 non-null     geometry
dtypes: float64(3), geometry(1), object(4)
memory usage: 1.6+ KB

Make this Notebook Trusted to load map: File -> Trust Notebook

6.1.8 Criar Mapa Temático em Folium

# Import packages
# import pandas as pd
# import geopandas as gpd
# import folium, folium.plugins 

# Mesmo Exemplo criar Plot dos Areas NUTS3 - Mostrando total de população 65+ no total de população
# Assegurar que foi efetuado Merge com Dados C2021
# GDF input: gdf_nuts3_2

# Selecionar Dados Portugal Continental: (copy() para criar nova coluna)
gdf_nuts3_sel = gdf_nuts3_2.loc[gdf_nuts3_2.NUTS1_x == '1'].copy()

# Novo Atributo
gdf_nuts3_sel['ind65'] = gdf_nuts3_sel.N_INDIVIDUOS_65_OU_MAIS/gdf_nuts3_sel.N_INDIVIDUOS

# Obter o centroid de toda a geometria
centroid = gdf_nuts3_sel.to_crs(epsg=4326).unary_union.centroid

# Criar Listagem com localização de latitude  longitude
center_map = [centroid.y, centroid.x]

# Criar Mapa e mostrar
folium_map = folium.Map(location=center_map, zoom_start=6, tiles='OpenStreetMap')

# Criar choropleth
choropleth = folium.Choropleth(
    geo_data=gdf_nuts3_sel,
    name='População 65+ mais',
    data=gdf_nuts3_sel,
    columns=['NUTS3', 'ind65'],
    key_on='feature.properties.NUTS3',
    fill_color='Reds', #
    fill_opacity=0.5,
    line_opacity=1.0,
    bins =8,
    legend_name='Relacao população 65+ no total de população'
)

folium.GeoJson(gdf_aces,
zoom_on_click = True).add_to(folium_map)

# Adicionar a mapa
choropleth.add_to(folium_map)

# Widget para controlar os layer visiveis            
folium.LayerControl().add_to(folium_map)

folium_map.save(r'data\xxmap.html')

folium_map

# Limpar Object da memoria
gdf_nuts3_sel = None